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by 

J.  A.  Tomlin 


Abstract 

This  report  gives  programmers  information  and 
documentation  for  LCPL— an  efficient  robust  program  for  solving 
Linear  Complementarity  Problems  by  Lemke's  Method. 
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DOCUMENTATION  OF  I.C1I 


T 


GENERAL 


Purpose 

LCP L is  a program  to  solve  linear  complementary  problems  (LCPs,' 
by  Letnke’s  method  l:i].  It  is  designed  to  handle  problems  with  a few 
hundred  rows  and  a sparse  coefficient  (M)  matrix.  Formally  an  LCP  is  a 
problem  of  the  form; 

Find  w and  z such  that 

w = p + Mz 

T 

w,  z > 0,  w z = 0 . 


Method 

The  code  uses  a product  form  implementation  of  Lemke's  method,  with 
an  efficient  inversion  routine  producing  a sparse  and  accurate  LU  decomposition 
of  the  basis  in  product  form.  The  basis  is  updated  by  the  standard  product 
form  method  with  accuracy  checks.  The  Harwell  scaling  subroutine  MCLsA 
using  Curtis  and  Reid's  scaling  method  is  incorporated.  A routine  for 
recovering  from  loss  of  feasibility  or  singularity  of  the  basis  is  provided. 

The  methodology  is  fully  described  in  reference  fl].  The  problem  is  stored 
in  the  form 

Iw  - Mz  - e£  = q 

where  e has  a +1  on  every  row  with  q < 0.  £ ("DUMMY  Z"  ' is  imme  iiately 

introduced  into  the  basis  to  give  a basic  feasible  almost  complementary 
solution.  The  complement  of  the  ejected  variable  is  then  introduced  and 
so  on.  The  algorithm  terminates  when  f ' ZTOLZE  (or  r is  non-basi?'  or 
when  a ray  is  encountered. 
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Linear  Complementarity 
Mathematical  Programming 
Quadratic  Programming 
Sparse  Matrices 

bounce  Language 

The  program  is  written  entirely  in  FORTRAN  IV  for  IBM  series  J6G  and 
570  computers.  It  is  WATFIV  compatible. 

Specification 

The  main  program  should  be  executed  as  a Job  step.  The  program 
input  is  described  in  the  "Users  Guide  to  LCPL"  [3],  All  input  occurs  in 
subroutines  INPUT  and  BASINP. 

Error  Indicators 

Error  indicators  and  other  diagnostics  and  messages  are  described  in  the 
s.ers  Guide  to  LCFL"  [3j.  They  will  also  be  indicated  in  the  documentation 
for  each  subroutine. 


S -b routines 


The  routines  making  up  LCPL  are 
BASINP 
BASOUT 
CHSOL 
CHUZR 
CLEAR 
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FT  RAN 

TTEROP 

INPUT 

INVERT 

LEMKE 

MCI  A (Harwell  Subroutine.. 

OUTPUT 

RECOVR 

SCALE 

SIIFTE 

SHIFT?. 

UNPACK 

UPBETA 

WRETA 

XCOL3 


| Program  Size 

The  total  length  of  the  program  is  1867  source  statements 
(including  comments  and  common  blocks,  etc.  . 

1 

«| 

Array  Sice 

The  ctandar i version  of  the  code  requires  a blank  common  array 
area  of  130,0'  bytes  f core.  This  is  of  course  dependent  on  the  setting 
of  maximum  i.roblem  j i mensions.  There  is  a labelled  common  requiring  (JO 

bytes.  This  is  not  variable. 
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Accuracy  and  Convergence 


The  code  has  been  specially  written  to  be  as  robust  as  possible 
(see  [4]).  If  the  matrix  M is  positive  semi-definite  or  has  positive 
principal  minors,  a recovery  procedure  is  provided  (RECOVR)  which  virtually 
assures  that  loss  of  feasibility  or  basis  singularity  can  be  overcome. 

This  facility  also  allows  the  user  to  start  such  problems  from  any  comple- 
mentary or  almost  complementary  basis.  Note  that  convergence  is  assumed 
when  the  dummy  variable  "DUMMi  !?'  falls  below  ZTOLZE.  It  need  not 
necessarily  leave  the  basis. 

Timing 

Using  the  FORTRAN  H compiler  with  OPT  =2  the  code  has  solved 
problems  of  order  100  in  about  h seconds  and  of  order  500  in  20  seconds 

on  an  IBM  370/168. 

References 

11  R.W.  Cottle  and  G.B.  Dantzig,  "Complementary  Pivot  Theory  of 

Mathematical  Programming,"  pp.  115-138  in  Mathematics  of  the  Decision 
Sciences  (G.B.  Dantzig  and  A.F.  Veinott,  Jr.,  Eds.},  Amer.  Math. 

Society,  Providence,  R. I.  (1968}. 

[2]  C.  E.  Lemke,  "Bimatrix  Equilibrium  Points  and  Mathematical  Programming," 
Management  Science  11,  pp.  681-689  (I9r5';. 

[3]  J.A.  Tomlin,  "Users  Guide  to  LCPL — A Program  for  Solving  Linear 
Complementarity  Problems  by  Lemke’ s Method,"  TR  SOL  7^-lf,  Systems 
Optimization  Laboratory,  Stanford,  Ca. , August  197f. 

[ 4]  J.A.  Tomlin.  "Robust  Implementation  of  Lemke’ s Method  for  the 
Linear  Complementarity  Problem."  TR  SOL  7P-2U,  Systems  Optimization 
Laboratory,  Stanford,  Ca. , September  1 <7(. 
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Further  Comment £ 


Since  virtually  all  communication  is  via  common  blocks,  and  all 
arrays  are  i n common,  the  contexts  of  both  labelled  and  blank  common 
are  described  in  the  next  section.  Furthermore  since  the  code  uses 
sparse  matrix  techniques  for  both  matrix  and  transformations,  and  these 


appear  in  many  subroutines,  these  techniques  are  also  described  in  detail. 


II.  COMMON  BLOCKS 


Labelled  Common 


COMMON /BLOCK/ 


This  common  block  contains  tolerances 


maximum  array  dimensions  and  literal  constants.  All  these  values  are 


set  in  BLOCK  DATA 


Tolerances  (REAL*4) 


ZTOLZE 


ZTOLPV 


Absolute  pivot  tolerance 


Absolute  tolerance  on  eta  elements 


ZTETA 


Relative  pivot  tolerance  used  in  CHUZR 


Dimensional  Parameters  (INTEGER*4) 


NRMAX 


Maximum  row  dimension 


NTMAX 


Maximum  number  of  eta  transformations 


NEMAX 


Maximum  number  of  eta  elements 


(4000) 


NAMAX 


Maximum  number  of  matrix  nonzeros 


These  values  must  be  changed  if  the  relevant  array  dimensions 


Literals  (INTEGER* 4 


All  literal  constants  begin  with  Q 


QA  is  'Abbb 


QH  is  ' Hbbb 


Cbbb 


QE  is  ' Ebbb 


Lbbb 


)0  is  'Obbb 


QU  is  ’Ubbb 


bbbb’  (blank),  QDUM1  is  ' DUMM' , QDUM?  is  'YbZb 


w 


Blank  Common 

Almost  all  communication  between  subroutines  is  by  means  of 
variables  and  arrays  in  blank  common.  (The  major  exception  is  the  Harwell 
scaling  routine  MC12A. ) We  indicate  the  dimension  with  each  array. 


\ 


REAIt'8  Arrays 

B (NHMAX)  The  given  right  hand  side 

X (ISSMAX)  Holds  the  values  of  the  basic  variables 

Y (NHMAX)  Work  region  (used  mainly  to  update  columns) 

YTEMP  (NHMAX ) Second  work  region 

The  k arrays  are  contiguous  and  considered  to  be  arrays  1 through  b by 
subroutine  3HIFTR  which  moves  information  between  them. 

E (NEMAX)  Contains  the  non-zero  eta  elements.  (They  are  not 

divided  by  the  pivot. ) 


REAL *8  Variables 

DSUM,  DPROD,  DY,  DE  and  DP  are  used  for  temporary  storage  in 
several  routines. 

ERMAX  contains  the  maximum  relative  error  found  in  the  right  hand 

side  during  CHSOL. 

J 

j 

REALvk  Array 

I A (NAMAX) 


Contains  the  non-zero  elements  of  the  matrix 
(I, -M, -e) 


hfc.  i ^ , 
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INTEGERS  Arrays 

ICNAM  (2*NRMAX*2,2  ) 


NAME  (20) 

NTEMP  (20) 

INTEGER*!*  Variables 
IBFRQ 

IFNEG 

IFSCAL 

INVFRQ 

I OUT 

IR0WP 

IR0WZ 

ITCH 

ITCHA 

ITCNT 


Contains  the  row  and  column  names. 
ICNAM(J,l)  and  ICNAM (,J, 2)  contain  an  8 
character  row  name  for  J < NR0W  and  a 
column  name  for  J > NROW.  The  " dummy- 
column"  name  "DUMMY  Z"  is  in  the 
2*NR0W+1  positions. 

Used  in  reading  input  cards.  The  first 
two  words  contain  the  problem  name  after 
INPUT. 

Temporary  storage. 


Iteration  frequency  for  saving  the  basis 
(if  KOUTB  specified) 

Switch  to  reverse  sign  of  matrix  elements. 
Standard  is  1 iyes).  Set  to  0 for  no. 
Switch  for  scaling.  Standard  is  0 (no). 
Iteration  frequency  for  inversion 
Output  option  switch.  Standard  is  0. 

Pivot  row 

Row  on  which  "DUMMY  Z"  pivots 
Iteration  frequency  for  calling  CHSOL 
Counter  of  iterations  since  last  call 

to  CHSOL 

Iteration  counter 


NRCW 


Number  of 


ows  in  the  problem 


NS  TNG 


NIJELEM 


NUETA 


INTEGERS  Arrays 


IA  (NAMAXi 


te(nemax) 


ISC (NRMAX+1 ,2 ) 


JH (NFMAX ! 


KINBAS  (2*NRMAX+2 ) 


LA  (2*NRMAX  - ' 


Number  of  singularities  found  by  INVERT 
Number  of  U eta  elements  in  INVERT 
Number  of  U etas  in  INVERT 


IA(l)  contains  the  row  index  for  the 
element  stored  in  All). 

IE  (1 1 contains  the  row  index  for  the  eta 
element  stored  in  Ell) 

ISC(I,1’  contains  an  integer,  say  R, 
which  defines  the  row  scaling  factor  lC1' 
for  row  I of  M. 

ISC (l ,2 ) contains  an  integer,  say  C, 
which  defines  the  column  scaling  factor 
16  for  column  I of  M. 

JHIl)  contains  the  number  of  the  vector 
Icolumn ) pivoting  on  row  I. 

KINBAS  1j)  is  zero  if  column  J is  nonbasic 
and  enual  to  the  row  on  which  it  pivots 
if  the  column  is  basic 

LA(l 1 points  to  the  first  nonzero  element 
of  column  I in  arrays  IAc  ) and  A1‘). 

The  last  nonzero  of  column  I is  pointed 
to  by  IA  (l-t-l ) -1. 
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J 


J 


LE (NTMAX ) 


LE(l)  points  to  the  firs*  nonzero  element 
of  transformation  I in  arrays  IE(*  ) and 


E(*  >.  This  is  always  the  pivot  row  and 
element  for  this  transformation.  The 
nonpivot  elements  follow.  The  last  non- 
zero of  transformation  I is  pointed  to 
by  LEd+l'-l. 


Changing  Array  Dimensions 

The  user  may  change  one  or  all  of  the  following: 
Maximum  row  size 

Maximum  number  of  nonzero  matrix  elements 
Maximum  number  of  eta  elements 
Maximum  number  of  eta  transformations. 


Row  Size 

The  value  NR MAX  in  labelled  common  must  be  changed  to  the  new 

maximum  row  size  in  BLOCK  DATA.  The  arrays  B(0,  X(*),  Y(*  ),  YTEMP(*  and 

JH(*1  must  be  dimensioned  to  the  new  value  of  NRMAX.  Similarly,  KINBASd) 

must  be  dimensioned  to  2*NRMAX+2,  as  should  LA(')  and  the  first  dimension 

of  ICNAM(.,2).  The  first  dimension  of  ISC(.,2)  should  be  NRMAX ->-1. 

In  addition,  the  equivalence  statement  in  subroutine  INVERT  must 

be  modified  so  that  HREd(l)  is  enuivalenced  to  YTEMP(—  NRMAX+l' . NRMAX  should 

2 

therefore  always  be  even.  Finally  BARRAY(- ) in  SHIFTR  (equivalenced  to  Bt,*'1'! 
should  be  dimensioned  to  4* NRMAX,  as  should  WSPACE(*)  in  SCALE  (equivalenced 

to  '/(')>. 
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Matrix  Elements 


t i._.  I 


ir^i 


The  value  LAi-IAX  in  labelled  common  must  be  changed  in  BLOCK  DATA. 
The  arrays  IA(')  . A • ) must  be  redimensioned  to  the  new  value  of  NAMAX. 


Eta  Elements 

The  value  of  'IE.MAX  in  labelled  common  must  be  changed  in  BLOCK  DATA. 
The  arrays  IE(*  EC*)  must  be  redimensioned  to  the  new  value  of  NEMAX. 


Eta  Transformations 

The  value  of  NTMAX  in  labelled  common  must  be  changed  in  BLOCK  DATA. 
The  dimension  of  LEi*  ! should  be  changed  to  NTMAX+2 . 


Overall  Arrav  Size 


Any  of  the  above  changes  will  change  the  length  of  blank  common. 
The  array  size  in  subroutine  CLEAR  should  be  changed  to  reflect  this. 


III.  INFORMATION  STRiiCTijRE." 


Matrix  Storage 

The  problem  matrix  (l,-M,-e)  is  stored  in  the  form  of  a threaded 
list  using  the  arrays  LA(*',  A(*  ) which  contain  respectively:  the 

pointer  to  the  first  entry  of  each  column  in  IA(-)  and  A(-^,  the  row 
indexes  of  the  nonzero  entries , and  the  values  of  the  nonzero.  All  elements 
of  a column  are  contiguous. 


Example. 

The  small  matrix 


1. 

-U. 

1. 

-3 

- 

1. 

2. 

■> 

i . 

would  be  stored  as  follows 

LA 

IA 

A 

1 

1 

1, 

2 

2 

1. 

X 

3 

1. 

k 

1 

-h. 

6 

3 

2. 

8 

2 

-3. 

10 

3 

1. 

12 

1 

-2. 

2 

-1. 

1 

-1. 

2 

-1. 

Note  that  there  is  one  more  entry  in  LAv‘)  than  there  are  columns. 
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I. 

I 


i 


This  points  to  the  position  after  the  last  entry  in  IA(* ! and  A(*  ) enabling 
this  last  entry  to  be  located.  To  unpack  say  column  J into  the  work 
region  Y(*  ) (previously  set  to  zero)  the  code  would  be  as  follows 


(Find  first  and 


KK  = LA(j+l)  - 1 


last  entry  of  column 


D0  100  I = LL,  KK 


(initiate  loop) 


IR  = IA(l) 


(Extract  row  index) 


y(ir^  = a (i : 


(Store  value  in  Y(-)) 


100  CONTINUE 


Transformation  Processing 


The  inverse  of  the  basis  is  maintained  as  a product  of  elementary 


transformations: 


where 


B =EtEt-l  "•  E2E1 


Each  h,  is  constructed  from  a vector  a as  follows: 


E,  to  obtain  y 


Hence  when  multiplying  a vector 


is  within  tolerance  of  zero 


Note  also 


that  it  is  only  necessary  to  store  the  a..  (not  the  p 


The  nontrival  columns  for  each  eta  are  stored  con' iguously  in  the 


same  way  as  the  columns  of  a matrix  using  the  three  • rrays 


Here  LE(-)  points  to 


■he  first  nonzero  of  each  transformation 


IE ( * ) contain  the  row  indices  of  the  nonzero  elements,  and  E."  1 contains 


the  nonzero  values  of  a 


The  pivot  of  each  transformation  is  identified  by  its  being  the 


first  element  in  IE(‘)  and  E('^  for  that  eta 


Thus  to  apply  transformation  K to  the  work  vector 


code  could  be  as  follow: 


KK  = I£(K*1)-1 
IPTV  - IE(LL) 

DY  = Y(IPIV) 

IF  (DABS  (DY)  .LT.  T0L)  G $)  T 0 1000 
DY  = Dif/EvLL) 

Y(IPIV)  - DY 

IF  (LL.GE.  KK  H 1 » 

LL  = LL+1 

D(Z>  500  I LL,  KK 

IR  = IE(l) 

Y(IR)  = Y(IR)  - E(l)*DY 
500  CONTINUE 


1000  CONTINUE 


Main  i rogram 

Calls  the  major  subroutines  of  the  code--INPUT,  .'PMKE  (the 
algorithm)  and  OUTPUT. 

The  action  of  the  main  program  is  determined  by  a parameter 
IACT  returned  by  IP PUT  as  follows: 

If  IACT  < 0 Terminate 
If  IACT  - 0 Call  LEMKE 
If  IACT  > 0 Call  OUTPUT 
OUTPUT  is  always  called  after  executing  LEMKE. 

After  OUTPUT'  the  program  calls  INPUT  again  to  look  for  a new  problem. 


BAS  IN  P 


Pends  a complementary,  or  almost  complementary  basis  from 


^■/  If  column  DUMMY  Z'  is  supplied  as  a column  its  accompanying 


row  name  is  checked  for  validity  and  the  column  replaced 
the  unit  vector  corresponding  to  the  row  name  in  the  basis. 

Input  Format 
NAME  Card 

Has  "NAME"  in  columns  1 - h and  the  basis  name  in  columns  15  - 22. 
If  this  name  does  not  agree  with  the  problem  name,  or  the  name  card 
is  missing,  a warning  is  given. 

Basis  Cards 

For  basic  columns  of  M the  column  names  are  given  (one  per 
card)  in  columns  5-12.  If  £ is  basic  there  must  be  a card  with 
"DUMMY  Z"  in  columns  5-12  and  a row  name  in  columns  15-22.  This 
row  name  is  the  name  of  the  w variable  £ replaces  in  the  basis. 

i 

i 

ENDATA  Card 

Has  "ENDATA"  in  columns  1-6. 

Note  that  if  the  basis  supplied  is  singular  or  infeasible 
the  program  will  enter  the  recovery  procedure  ("RECOVR").  If  naming 
errors  are  encountered  in  reading  the  basis  the  program  will  start  with 
the  partial  basis  found  (which  will  almost  certainly  trigger  a call  to 
RECOVR ) . 
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STARTING  BASIS  WILL  BE  READ  PROM  PILE  'number'. 


■ 

I 

, 


Message  that  KOUTB  was  specified. 


WARNING-BASIS  NAME  CARD  MISSING 


WARNING-BASIS  HAS  NAME  ' name' -DIFFERENT  FROM  PROBLEM 

Warning  messages  on  reading  a basis.  May  not  be  the  correct 
basis  for  the  problem.  Run  continues. 

NO  MATCH  FOR  VECTOR 

The  basis  file  has  a vector  name  not  found  in  the  problem. 
Vector  ignored  and  run  continues.  Will  almost  certainly  yield  singular 
or  infeasible  basis  and  a call  to  the  recovery  routine. 

NO  MATCH  FOR  ROW  'rownam'  TO  SWAP  WITH  DUMMY  Z 


j The  row  name  'rownam'  cannot  be  found  when  reading  the  basis 

, and  the  "DUMMY  Z"  variable  is  not  put  in  the  basis. 

The  run  will  continue. 

SPURIOUS  CARD  IN  BASIS  FILE-' text' -ABORT 

i 

, An  unrecognizable  card  in  the  basis  file  - fatal  error. 


!>i 


f 


BASOUT 


Writes  the  names  of  the  basic  structural  (m)  columns,  as  well 
as  " DUMMY  Z"  and  the  row  on  which  it  pivots  to  file  KOUTB. 


Formal  Declaration 


SUBROUTINE  BASOUT 


Method 


0/  If  KOUTB  not  specified  exit.  If  KOUTB  > ? rewind. 

1/  Write  NAME  card  (see  BASINP  for  format) 

2/  Write  names  of  basic  structural  columns 

3/  If  "DUMMY  Z"  basic  check  that  NOIJC  defined.  If  so  find  unit 
vector  "DUMMY  Z"  replaces  and  write  their  names. 


Messages 

NONCOMPLEMENTARY  PAIR  NOT  DEFINED,  BUT  DUMMY  Z IN  BASIS,  UNRESOLVABLE 
ERROR -NO  BASIS  SAVED. 

The  solution  is  supposedly  complementary  but  DUMMY  Z is  basic. 
Highly  unlikely  unless  program  corrupted. 


BASIS  WRITTEN  ON  FILE  'number'  at  ITERATION'  number’ 
The  basis  has  been  written  on  file  KOUTB. 


CHUZR 


Selects  the  variable  to  leave  the  basis  and  the  pivot  row  for 


the  incoming  variable. 


Formal  Declaration 


SUBROUTINE  CHUZR  (iFPV) 


Method 


Parameter:  IFPV  (Output)  set  to  -1  if  infeasibility  observed 

set  to  0 if  unacceptable  pivot  found 

set  to  1 if  acceptable  pivot  found. 


Performs  a ratio  test  on  the  elements  of  the  updated  increasing 


vector  Y(* ) and  the  solution  X(‘ ) to  determine  the  pivot  row,  by 
P.M.J.  Harris'  two  pass  technique.*  The  solution  X i * ) is  also  examined 
for  in feasibility.  The  technique  is  as  follows; 


_ . X(l)  * ZTOLZE 

Find  DPSI  = min ; — : 

Y(l ) 


for  Y(l)  > ZTOLPV 


If  X(l)  < -ZTOLZE  set  IFF/  = -1  and  exit 


Find  max  Y(l ! 

DTHETA < DPSI 
Y ( I ) > ZTOLPV 


x(i ) 

where  DTHETA  = rrrn' 


and  choose  row  IROWP  which  gives  maximum. 

5/  If  Y(lROWP'  > ZTOLRP*max I Y(l ) | set  IFPV  = 1,  otherwise  set 

IFPV  = 0. 

v/  If  X(lROWP)  <0  set  to  zero. 

Reference:  P.M.J.  Harris,  "Pivot  Selection  Methods  of  the  Devex  LP  Code," 
Mathematical  Programming  Studies  4,  (1975'..  .30-5?. 
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Message 


FEAS  LOST  ON  ROW  'row  no."  VAR  'variable  name'  'value' 

A loss  of  feasibility  is  detected  when  attempting  to  choose  a 
pivot  row.  Call  to  RECOVR  is  triggered  in  LEMKE. 

CLEAR 

Sets  all  of  blank  common  to  zero. 

Specification 

SUBROUTINE  CLEAR 


Method 

Blank  common  is  defined  as  a single  array  - ARRAY(*)  - which  is 
set  to  zero  in  a loop.  The  length  of  ARRAY (*  ' should  be  the  number  of 
single  words  in  blank  common  in  all  other  routines. 

FT  RAN 

p'/'f-s  TL‘  wot  ' vector  Y(*  by  the  current  inverse  representation 
eta  file  . 

Formal  Declaration 

SUBROUTINE  FTRAN  (iPAR) 

Parameter:  TRAP,  (input1  is  set  to  0 if  Y * is  to  Ve  updated 

by  the  entire  eta  ri le  (eras  1 through  NETA ' . IPAR  is  set 
to  1 to  avoid  updating  '{('  ) by  above  bump  etas  in  INVERT 
(Y(')  is  updated  by  etas  NUET.V  1 through  NETA. 


Method 


i 


The  standard  LP  FTRAN  procedure  is  followed.  Phe  data  structure 
of  the  eta  file  is  described  in  the  section  "Transformation  Processing." 

INPUT’ 

Reads  the  non-standard  parameter  settings  for  the  run  (if  any' 
and  the  problem  matrix  M and  right-har.d  side  q in  a slightly  modified 
MPS  format. 

Specification  and  Parameters 

SUBROUTINE  INPUT  (TACT) 

Parameter:  IACT  (Output'  Set  to  -1  if  the  subroutine  detects  the 

end  of  the  input  stream,  or  if  a fatal  input  error  is  dis- 
covered. Set  to  0 if  the  problem  is  properly  posed  and 
does  not  have  a trivial  solution.  Set  to  1 if  problem  read 
has  a trivial  solution  (i.e.  unit  or  given  basis  is  feasible 
and  complementary ' . 
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Program  Parameter  Input 

Parameters  are  changed  from  their  default  value  for  each 
problem  by  means  of  a FORTRAN  "NAMELIST"  input  from  the  card  reader. 
The  first  card  must  start  with  &PARAM  in  column  2 or  later*  each 
card  must  begin  with  a blank, and  the  list  must  be  terminated  by  &END. 


ol 


I; 

r 

| 

I 
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Examples. 


cols 


&PARAM  NQUAD  61+,  IF3CA1>1,  EiiD 

&PARAM  4END 


The  first  example  specifies  that  the  problem  is  a quadratic  program 
in  the  first  6k  variables  and  calls  for  scaling  of  M.  The  second 
example  leaves  all  parameters  at  their  default  values. 

On  normal  completion  of  input  we  return  to  the  main  program  with 
the  packed  problem  matrix  [l,-M,-e]  stored  in  LA(*),  IA(*  A \ • ' and 
the  right-hand  side  in  B(*  ).  All  the  row  and  column  names  are  held  in 
ICNAM(.,.).  The  initial  basis  is  defined  by  JK(*  ' and  KINBAS  • . The 
row  on  which  "DUMMY  Z"  has  pivoted  is  stored  in  IROWZ  and  the  variable 
thus  eliminated  from  the  basis  is  stored  in  NONC. 
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We  classify  the  parameters  into  the  following  groups: 
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Hu n Control; 

ITRLIM  (default  = 99999) 

Maximum  number  of  iterations  for  this  run  of  this  problem. 

INVFRQ  (default  = 30 ) 

Inversion  frequency 

ITCH  (default  = INVFRQ) 

Iteration  frequency  for  checking  the  relative  error  in  the 
current  basic  solution.  Defaulted  to  do  so  just  before  normal 
inversion.  (INVERT  automatically  carries  out  this  check 
after  inversion. ) 

IBFRQ  (default  = 99999) 

Iteration  frequency  for  writing  the  current  basis  to  the 
file  specified  by  KOUTB  (not  done  unless  KOUTB  is  specified 
to  be  non-zero).  If  KOUTB  is  properly  specified  the 
final  basis  will  be  saved  regardless  of  the  value  of  IBFRQ. 


Input : 

KINP  (default  = 5) 

FORTRAN  logical  unit  number  for  reading  the  problem. 
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IFNEG 


(default  = l) 


Determines  whether  to  change  the  sign  of  the  input  matrix  M. 

With  the  default  setting  the  user  supplies  the  matrix  M and 
-M  is  stored  to  collect  all  variables  on  the  left-hand  side 
as  in  (6).  If  the  user  supplies  -M  then  set  IFNEG  = 0. 

IFSCAL  (default  = 0) 

This  parameter  is  set  to  1 if  the  matrix  M is  to  be  scaled. 

If  scaling  is  performed  the  solution  (and  quadratic  objective 
function  value)  are  descaled  to  give  their  true  values.  Usually 
scaling  will  not  be  required. 

IFALLE  (default  = 0) 

If  the  parameter  is  set  to  1 the  dummy  column  e will  be 
constructed  with  a 1 in  every  position.  Not  recommended  for 
normal  use  since  this  increases  the  density  of  the  transformations. 

KINB  (default  = 0) 

FORTRAN  logical  unit  number  for  reading  a basis.  With  the 
default  value  of  0 no  basis  reading  will  be  attempted.  If 
KINB  is  greater  than  7 the  unit  is  rewound  before  reading. 

NQUAD  (default  = o) 

If  NQUAD  is  non-zero  the  problem  is  assumed  to  be  a quadratic  program 
set  up  as  in  the  "Users  Guide".  In  this  case  only  the  NQUAD 
columns  of  M corresponding  to  the  quadratic  primal  variables 
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need  be  supplied;  the  -A  column  will  be  constructed  from 


the  last  NROW-NQUAD  rows  of  the  column  supplied.  If  M is  non 


square  and  NQUAD  is  no' 


correctly  specified  a fatal  error 


(default  = 0) 


With  the  default  setting  all  the  w and  z variables,  name 


and  values, are  printed.  If  I0UT  is  set  to  a non-zero  value 
only  the  basic  variables  are  printed  (in  sorted  order)  together 


with  the  row  names  and  the  original  right  hand  side  q 


(default 


NQUAD 


If  NQUAD  has  some  non-zero  value  n the  first  n rows 


and  columns  of  M are  assumed  to  be  the  Hessian  of  a quadratic 


form.  The  value  of  the  quadratic  form  will  be 


printed  out  after  the  w and  z variable  values  at  solution 


time.  (See  also  Input  above. ) 


(default  = 0) 


KOUTB 


FORTRAN  logical  unit  number  for  writing  the  basis  every  IBFRQ 


is  exit  put  if  KOUTB 


iterations  and  on  termination.  No  basi 


is  zero.  Unit  KOUTB  is  rewound  before  and  after  writing 


if  KOUTB  is 


Tolerances: 

ZTOLZE  (default  = lO"** ) 

Zero  tolerance  for  testing  feasibility.  Basic  variables  are 
non-negative  if  they  are  > -ZTOLZE.  Complementarity  is 
declared  if  t,  ("DUMMY  2')  <_  ZTOLZE. 

ZTOLPV  (default  = 10'6) 

Absolute  pivot  tolerance. 

-8 

ZTOLRP  (default  = 10  ) 

Relative  pivot  tolerance  while  iterating.  A pivot  is  accepted 
if  its  absolute  value  is  greater  than  the  largest  absolute 
value  in  the  updated  column  multiplied  by  ZTOLRP.  If  this 
test  fails  INVERT  is  called  and  the  same  variable  is  considered 
again.  If  the  chosen  pivot  element  still  fails  the  test  it 
is  accepted.  Note  that  if  ZTOLRP  is  large  this  could  lead 
to  inverting  at  every  iteration.  If  ZTOLRP  is  too  small  the 
test  becomes  ineffective. 

ZTETA  (default  = 10'1?) 

Tolerance  on  the  absolute  value  of  transformation  elements. 
Larger  values  may  give  sparser  transformations  but  at  great 
cost  to  accuracy  and  stability. 

•7 

ZTOLDA  (default  = 10  ) 

Tolerance  on  the  absolute  value  of  entries  in  the  matrix  M. 

Any  value  smaller  than  this  is  disregarded. 


Problem  input 


The  M and  q defining  be  problem  are  read  from  unit  KINP 
in  slightly  modified  "MP3  format..''  Phe  card  images  required  are: 

NAME  Card 

This  has  "NAME"  in  columns  1 to  1.  and  the  (up  to  H characters 
problem  name  In  columns  15-22.  'ibis  card  is  optional  bu  . highly 
desirable  (particularly  for  checkin'.  against  basis  names'. 

ROWS  jf_ard 

Has  "BOWS"  in  columr:s  1-1. 


Row  dames 

Each  row  is  assigned  a (unique)  name  of  up  to  8 characters,  one 
per  card,  in  columns  5 -12.  'Embedded  blanks  are  allowed  (unlike  MPS). 
b'ote  that  there  are  no  row  * ypes  as  in  MK'  and  tha  any  supplied  will 
be  ignored.  Column  1 must  be  blank. 

Cf‘LU’W;  :ard 

Has  ' COLUMTin"  in  columns  1-7. 

Matrix  Element 

The  non-zero  elements  of  M are  supplied  by  column.  All  the 
f lements  of  a column  must  be  together.  Each  eol'anr,  is  assigned  a 
'.unique)  name  of  up  to  eight  characters.  The  forma  is: 
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cols 


15  - 22 


1-4  5-12 

blank  column  row 

name  name 


25  - 36  40  - 47  50  - 6l 


element 

value 


second 

row 

name 

(optional) 


second 

element 

value 

(optional) 


Again  embedded  blanks  are  allowed  in  column  names. 


RHS  Card 

Has  "RHS"  in  columns  1 - 3> 

Right  Hand  Side  Elements 

The  right  hand  side  vector  (q)  may  be  given  a name,  which 
should  be  different  from  any  row  or  column  name.  The  elements  are  given 
in  the  same  format  as  the  matrix  elements.  Note  that  only  one  right 
hand  side  may  be  supplied.  If  an  attempt  is  made  to  input  more  than 
one,  the  non-zeros  from  later  right  hand  sides  will  over-write  the 
earlier  values. 


SltDATA  Card 

Has  "ERRATA"  in  columns  1-6. 

Sample  input  is  shown  in  section  4.1. 

It  is  important  to  note  that  the  "w"  variables  are  assigned 
the  row  names  of  M,  and  the  "z"  variables  the  column  names.  The  dummy 
column  e is  assigned  the  name  " DUMMY  Z",  as  is  the  dummy  variable  t. 
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Messages 


SPLIT  VECTOR  'name' 

A column  with  the  same  ’name'  was  encountered  earlier.  The  input 
deck  is  probably  corrupted.  Input  continues  however 

NO  MATCH  FOR  ROW  ' rovname ' AT  COLUMN  'colname' 

Column  'colname'  contains  a nonzero  on  an  unidentifiable  row 
'rowname'.  This  is  a fatal  error. 


NUMBER  OF  MATRIX  ELEMENTS  'number'  EXCEEDS  NAMAX -ABORT 

The  program  dimensions  for  A^)  and  LA  ’ ) must  be  increased 
to  solve  the  problem.  Fatal  error.  Standard  NAMAX  is  LOGO. 

MATRIX  NOT  SQUARE  AND  NQUAD  INCORRECTLY  SPECIFIED- PROBLEM  ABANDONED. 

The  number  of  structural  columns  supplied  is  not  equal  to 
either  NF.OW  or  the  specified  value  of  NQUAD.  M cannot  therefore  be 
properly  constructed  to  be  square.  Fatal  error. 

NQUAD  COLUMNS  ONLY  SUPPLIED- -WIT, L ATTEMPT  TC  FORM  ’number'  EXTRA  COLUMNS 
OF  CONSTRAINT  MATRIX  TRANSPOSE 

The  problem  is  assumed  to  be  a quadratic  program 
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DUMMY  COLUMN  WITH  'number'  ELEMENTS  CONSTRUCTED 

There  are  negative  elements  in  the  q vector  and  a dummy  column 
has  been  constructed  for  "DUMMY  Z". 


PROBLEM  HAS  TRIVIAL  SOLUTION 

Either  q is  non-negative  or  the  given  starting  basis  yielc 
a feasible  complementary  solution.  No  dummy  column  is  constructed 
and  the  solution  is  printed  out. 


EASIS  GIVEN  NOT  FEASIBLE,  WILL  ENTER  RECOVP 

The  given  basis  (complementary  or  almost  complementary)  is 
infeasible.  The  recovery  routine  is  called  to  construct  a new  basis 
and  dummy  column. 


Subroutines  Called 


CLEAR 


BAS  IN  P 


FFCOVR 


INVERT 


XCOLS 


INVE'-  . 


Performs  an  LU  decomposition  of  the  basis,  using  a ""erit 


Method"  to  preserve  sparsity,  and  writes  the  product  form  of  U~  n" 


to  the  eta  file. 


7 / 


Formal  Declaration 


SUBROUTINE  INVERT 


* 

Method 

The  basis  is  (notionally)  permuted  to  the  form 


' L 

’ L1 

A 3 

and  then 

C U E 

- C E L2  - 

. A B 

where  is  lower  triangular  (the  above  the  bump  columns)  and  U^ 

upper  triangular,  permuted  from  lower  triangular  Ln  (the  below 

cl 

bump  columns).  B is  the  bump. 

The  resulting  factorization  is  of  the  form 


where  B has  triangular  factors  L^U^. 

The  routine  proceeds  by  pulling  out  all  the  above  and  below 
bump  etas  and  writing  them  out,  then  writing  the  bump  columns  in 
"merit"  order.  The  bump  is  then  decomposed.  The  resulting  L and 
U etas  are  merged  by  SIIFTE. 

* 

J.  A.  Tomlin,  "Pivoting  for  Size  and  Sparsity  In  L.P.  Inversion  Routines," 
J.  Inst.  Maths.  Applies.  10,  pp.  289-295  (1972). 
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On  completion  the  right-hand  side  is  placed  in  '/>•),  updated 
and  put  in  X(‘ )•  CHSCL  is  called  to  check  relative  error  in  X(* 
if  this  is  too  large  reinversion  is  performed  with  a larger  relative 
pivot  tolerance. 

Since  INVERT  is  the  most  intricate  routine  we  describe  it  more 
fully  as  follows: 


Local  Arrays 


H-region  Contains  sequence  numbers  of  vectors  already  pivoted  on. 

(HREG ) Contains  -1  - Row  Count  for  unas signed  rows 

Contains  0 for  rows  below  the  bump  and  not  yet  pivoted  in. 


V- region  Contains  sequence  numbers  of  columns  not  yet  pivoted  into  rows. 

(VSEG)  It  is  divided  into  k parts  ( some  of  which  may  be  vacuous) 

as  follows : 


Part  1 contains  sequence  numbers  of  columns  in  the  bump. 
Part  2 is  spare. 

Part  3 contains  sequence  numbers  of  structural  columns 
pivoting  below  the  bump 


Part  4 


contains  sequence  numbers  of  slacks  and  artificials 
pivoting  below  the  bump 


M- region  Contains  information  about  those  rows  mentioned  in  the  V-region 

(MREG)  It  is  divided  into  4 parts  exactly  like  the  V-region  as  follows 

Part  1 contains  the  merit  measure  of  columns  in  the  bump 

Part  2 is  spare. 
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Part  3 contains  the  pivot  rows  of  vectors  mentioned  in 
Part  3 of  V-region. 

Part  4 contains  the  pivot  rows  of  vectors  mentioned  in 
Part  4 of  V-region 

Note  that  MREG(- ) is  equivalenced  to  YTEMP ( • ' , VREG ( • to  Xt-  ) arid 
HEEG(*  ) to  YTEMP  r NRMAX+1  > . Thus  the  last  equivalence  must  be 
changed  if  NR  MAX  is  changed.  The  three  arrays  are  INTEGER *2  and 
of  dimension  NRMAX. 

OUTLINE  FLOW  CHART  FOR  INVERT 

1.  Scan  H-region 

If  a row  contains  a slack  or  artificial, 

put  sequence  number  into  Part  4 of  M-region  and  V-region 
Otherwise  put  sequence  number  into  Part  1 of  V-region. 

In  either  case  replace  H-region  by  -1. 


2.  Take  vectors  in  Part  4 of  M-region 

Zero  out  those  entries  in  the  H-region  mentioned  in  Part  4 

of  M-region,  and  pivot  these  (unit)  vectors  on  their  own  rows. 


Scan  Part  1 of  V-region 

Establish  row  counts:  Search  for  entries  in  rows  where  the  H-value 

is  negative  and  subtract  1 from  this  H-value.  If  there  are  no 
such  rows,  the  vector  cannot  be  introduced.  Write  message  to 
declare  MATRIX  SINGULAR  arid  assign  complementary  unit  vector  to 
basis.  Replace  vector  by  the  bottom  entry  in  Part  1 of  V-region. 
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If  there  is  only  one  such  low,  the  vector  will  pivot  below  the 
bump  in  this  row,  if  the  pivot  is  larger  than  V LPV. 

replace  vector  by  the  bottom  entry  in  Par!  1 of  V-region 
put  sequence  number  into  Part  3 of  V-region 
put  pivot  into  Part  3 of  M-region 
put  0 into  H-region 
Otherwise  go  on  to  next  vector. 


4.  Repeat  following  step  until  no  more  vectors  can  be  removed,  from  bump. 
Scan  Part  1 of  V-region 

If  a vector  has  an  element  in  a row  where  the  H-value  is  -2  then 
this  is  the  only  column  with  an  element  in  that  row  and  the 
vector  will  pivot  above  the  bump. 

Write  ETA  to  pivot  in  this  row. 

Add  1 to  negative  entries  in  rows  where  this  vector  has 
a coefficient. 

Put  sequence  number  into  --region. 

Replace  vector  by  bottom  entry  in  Pari  1 f -region. 
Otherwise  test  the  count  of  entries  in  any  rows  where  the  H-value 
is  negative. 

If  there  are  no  such  rows,  the  vector  cannot  be  introduced.  Write 
message  to  declare  MATRIX  SINOtJIAB  and  assign  complementary  unit 
vector  to  the  basis.  Replace  vector  by  the  bottom  entry  in 
Part  1 of  V-region. 
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If  there  is  only  one  such  row,  the  vector  will  pivot  below  the  bmp 
in  this  row,  if  the  pivot  is  larger  than  ZTOLiV 

Replace  vector  by  the  bottom  entry  in  Part  1 of  V-region. 

Put  sequence  number  into  Part  3 of  V-region. 

Put  pivot  into  Part  3 of  M-region. 

Put  0 into  H-region. 

If  there  is  more  than  one  such  row,  calculate  the  Merit  to  be  put 
into  the  M-region  as  the  sum  of  (row  count  -1J  over  all  rows  with 

hreg(i)  < 0. 

Sort  rows  in  Part  1 of  V-region  into  acsending  merit  order. 

Write  out  below  bump  etas  for  nonunit  vectors. 

Take  vectors  in  Part  1 of  V-region  in  turn. 

Unpack  entries  in  vector  into  Y i • ) 

Co  FTRAN  on  vector. 

Set  pivot  tolerance  ZTOLY  to  ZTOLPV 

Select  pivotal  row  as  the  one  with  coefficient  above  ZTOLY  in 
magnitude  that  has  the  smallest  row  count.  If  selected  pivot  has 
absolute  value  less  than  ZTRELy  (.maximum  eligible  Y 1 1 ) - reset 
ZTOLY  to  this  value  and  repeat  row  selection. 
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Write  L and  U etas  pivoting  on  this  row. 

Increase  row  counts  in  all  rows  where  ETA  nonzero  value  and 


H-region  is  negative  by  (pivotal  row  count  -2),  i.e. 
H-region  by  (pivotal  H-value  + 3). 


increase 


If  no  rows  have  a coefficient  above  E'fOLPV  in  magnitude  then  this 
vector  cannot  be  introduced.  Write  message  to  declare  MATRIX 
SINGULAR  and  assign  complementary  unit  vector  to  basis 

?.  Merge  L and  U etas  (calling  SHFTE) 

8.  If  singularities  were  detected,  print  message  and  exit 
Su  Update  solution,  i.e.  compute  X ( •) 

10.  Print  statistics 

11.  Check  relative  error  in  solution  (using  CHSOL)  and  if  the  value  is 
excessive  repeat  INVERT  with  ZTREL  increased  to  0-5- 

Messages 

MATRIX  SINGULAR -VECTOR  'column  name'  REMOVED 

INVERT  has  detected  a singularity,  bee  below. 

INVERT  FAILS.  WILL  ATTEMPT  TO  RECOVR 

One  or  more  singularities  were  detected  by  INVERT.  A call  to 
RECOVR  is  triggered. 

INSUFFICIENT  ETA  SPACE  TO  INVERT-RUN  A BAN LONER 
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INCREASE  ETA  SPACE  AND  NEMAX 

The  arrays  IE l'  ) and  E(* ) in  the  program  must  be  increased  from 
their  current  dimension  (standard  value  is  NEMAX  = 8000). 

Subroutines  Called 
WRETA 
FTRAN 
SHIFTR 
SHFTE 
CHSOL 


ITEROP 

Prints  the  iteration  log  (see  "Users  Guide  to  LCPL"). 


Formal  Declaration 

SUBROUTINE  ITEROP  (IPAR) 

Parameter:  IPAR (Input).  When  IPAR=0  print  headings  for  iteration  log. 

When  IPAR  j-  0 print  a line  of  the  iteration  log. 

Information  Printed 

ITCOUNT  = iteration 

A almost  complementary 
STATUS  = C complementary 

U unbounded  almost  complementary  ray 

DUMMY  Z = value  of  the  "dummy  variable" 

VECIN  name  of  column  entering  basis 
VECOUT  - name  of  column  leaving  basis 
NETA  - number  of  etas  in  the  eta  file 

NELEM  --  number  of  elements  in  the  eta  file. 
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Formal  Declaration 

SUBROUTINE  LEMKE 

Method 

0/  It  is  assumed  that  LEMKE  is  called  with  a feasible,  almost 
complementary  basis  found  in  INPUT. 

l/  INVERT  is  called.  If  the  basis  is  singular  a call  to  RECOVR 

is  triggered.  When  (or  if)  a new  nonsingular  almost  complementary 
basis  has  been  found  in  RECOVR.  the  INVERT  routine  is  called 
again. 

2/  After  INVERT  headings  are  printed  out  for  the  iteration  log 
(by  a call  to  ITEROP(o)). 

3/  The  row  IROWZ  on  which  "DUMMY  Z"  pivot  is  found. 

b/  Test  the  variable  NONC  which  has  just  left  the  basis  to  see 
if  it  is  a w variable  (NONC  < NROW),  in  which  case 
JCOLP  = NROW  + NONC  is  to  enter  the  basis,  or  a z variable 
(NONC  > NROW/  in  which  case  JCOLP  = NONC  - NROW  is  to  enter 
basis . 

5/  Column  JCOLP  is  UNPACKed  and  FTRANed. 
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6/  CHUZR  (iFPV)is  called.  If  IFFV=1  proceed  to  step  7/. 

If  IFPV  < 0 go  to  step  l/and  INVERT  unless  inversion  has 
just  been  performed.  If  ITSINV=0  and  IFPV--1  the  solution 
is  infeasible--go  to  RECOVR  and  step  1/.  If  ITSINV=0  and 
IFFV=0  the  pivot  row  IROWP  chosen  in  CHUZR  must  be  accepted. 

7/  Bring  column  JCOLP  into  the  basis,  pivoting  on  row  IROWP 
and  ejecting  the  column  JH (IROWP).  Redefine  NONC. 

8/  Update  X(*)  and  print  a line  of  the  iteration  log. 

9/  If  ’’DUMMY  Z"  (column  NCOL'  has  left  the  basis  or  its  value 
is  less  than  ZTOLZE  go  to  13/  and  declare  a complementary- 


solution 

10 / Write  a new  eta 

11/  If  it  is  time  to  do,  write  out  the  basis,  call  CHSOL  or 
call  INVERT  (i.e.  go  to  1 /). 

12/  If  ITCNT  is  > ITRLIM  terminate  and  go  to  13/ 

13/  Write  the  final  line  of  the  iteration  log  giving  terminal 
status,  the  final  basis,  and  return  to  the  main  program. 


Message 

ITERATION  LIMIT  EXCEEDED 


Subroutines  Called 

RECOVR 

UNPACK 

UPBETA 

BAS OUT 

INVERT 

FTRAN 

WRETA 

ITEROP 

CHUZR 

CHSOL 

1*2 


MC12A 


Harwell  subroutine  for  scaling  a square  matrix  using  the  method 


of  Curtis  and  Reid. 


Formal  Declaration 

SUBROUTINE  MC12A  (A,  IND,  IP,  N,  NP,  DIAG,  RES,  IS) 


The  parameters  and  the  methods  are  fully  explained  and  described  in 

A.R.  Curtis  and  J.K.  Reid,  "Fortran  Subroutines  for  the  Solution 
of  Sparse  Sets  of  Linear  Equations,"  Report  AERE-R  6841+ , 
Theoretical  Physics  Division,  AERE,  Harwell,  England  (l97l). 


Note 

Very  minor  changes  have  been  made  so  that  nested  DO-loops 
do  not  end  on  the  same  statement,  thus  achieving  WATFIV  com- 
patibility. 


OUTPUT 

Write  the  solution  found  by  LEMKE  to  FORTRAN  logical  unit  KOUTP. 
Formal  Declaration 

SUBROUTINE  OUTPUT 

Method 

0/  Write  the  problem  name  at  the  top  of  a new  page, 
l/  If  I0UT=O  go  to  6/  (Standard 
2/  Write  a heading 

3/  Sort  the  basic  variable  values  X(-  ) into  matrix  order  using 
the  sequence  numbers  in  JH(') 

1+3 





B - — v. 


4/  For  1=1,  NROW  write  the  name  and  value  of  a basic  variable 
(in  sorted  order),  the  row  name  and  the  value  of  the  original 
right-hand  side  element  B(l). 

5/  Go  to  9/. 

6/  Write  a heading 

7/  For  1=1,  NROW  write  the  name  and  value  of  the  w.  and  z 

x 1 

variables 

8/  If  "DUMMY  Z"  is  basic  write  its  name  and  value  under  the 
Z variables. 

9/  Compute  and  write  the  value  of  the  quadratic  objective  value 
if  NQUAD  has  been  specified. 


Note  If  IFSCAL=1,  all  the  solution  values  must  be  descaled  before 

pivoting  and  before  computing  the  quadratic  objective  function 
(see  SCALE  for  details  of  scaling  factors). 


REOOVR 


Construct  a new  nonsingular  almost  complementary  feasible 


basis  (if  this  is  possible)  from  the  old  basis  using  as  many  of  its 


structured  columns  as  nossible. 


Formal  Declaration 


SUBROUTINE  RECOVR 


44 


Method 


1 / Femove  "DUMMY  Z"  from  the  basis  and  replace  it  by  NONC  to 
give  a complementary  basis 

2/  Check  for  correct  number  of  basic  vectors  and  reset  JR(*  ) 
to  include  them.  Count  number  of  structural  columns 
(columns  of  M)  in  basis,  denoted  NSTB. 

3/  Attempt  to  invert  basis.  If  singularity  was  found  repeat 
until  for  a maximum  of  NSTB  tries  until  a nonsingular  basis 
is  found  (INVERT  automatically  ejects  dependent  structural 
columns  and  replaces  them  by  their  complementary  unit  column) 

h/  Examine  the  basic  solution  in  array  X(').  Whenever  a negative 
X(l)  is  found  subtract  the  vector  IV=JH(l)  from  region 
Y(« ) (initialized  to  zero).  Record  row  IROWZ  containing 
the  most  negative  entry 

5/  If  solution  is  nonnegative  (i.e.  IR0WZ=0)  declare  termination 
since  the  solution  has  been  constructed  to  be  complementary. 

6/  If  IROWZ  is  nonzero,  pack  the  column  in  LAI* ) and  A(*  ). 

Write  a message  that  this  has  been  done. 

7/  Define  NONC  as  JH(IROWZ)  and  change  TH(*  ) and  KINBAS v * ) so 
that  the  new  basis  has  "DUMMY  Z"  (i.e.,  column  NCOL) 
basis  in  place  of  NONC.  Return.  (INVERT  will  be  called 
again  on  return  for  the  new  almost  complementary  feasible 
basis . ) 


Ur> 


Messages 


FATAL  ERROR-’ number'  BASIC  VECTORS 

The  problem  does  not  have  NROW  basic  vectors.  Probably  due 
to  corrupted  program  or  split  vectors. 

COMPLEMENTARY  BASIS  FOUND 

RECOVR  has  found  a new  complementary  basis. 

FEAS  COMP  BASIS  ARRIVED  AT  IN  RECOVR 

RECOVR  has  (fortuitously)  found  the  complementary  basis  it 
has  constructed  to  be  feasible.  Problem  solved. 

NEW  DUMMY  COLUMN  CONSTRUCTED 

RECOVR  has  constructed  a new  dummy  column  and  determined  a 
pivot  row  which  will  give  a feasible  almost  complementary  solution. 
Returns  to  normal  algorithm. 

Subroutines  Called 
INVERT 


SCALE 

Sets  up  arrays  and  parameters  for  calling  the  Harwell  scaling 
subroutine  MCI? A and  then  applies  the  scaling  factors  returned  by 
MC12A  to  the  matrix  (m)  and  the  right-hand  side  stored  in  B(*). 


kC 


|,a 

t 


j 


Formal  Declaration 


SUBROUTINE  SCALE 

Special  Techniques 

Scaling  factors  are  applied  by  adding  integers  to  the  float- 
ing point  exponents  of  A(* ) (single  precision)  and  B(')  (double  pre- 
cision). This  is  accomplished  by  means  of  equivalences.  Single  and 
a double  precision  variables  AVAL  and  BVAL  are  defined  which  are 
equivalenced  to  L0GICAL*1  variables  IU  and  IV.  IU  and  IV  are  thus 
the  exponents  of  AVAL  and  BVAL  with  the  standard  bias  of  64  (decimal) 
added.  Any  matrix  element  A(l),  or  right-hand  element  B(l),  is  thus 
stored  in  AVAL  or  BVAL  and  the  exponent  is  extracted  or  modified  by 
examining  or  changing  IU  or  IV.  The  modified  AVAL  or  BVAL  then 
overwrites  A(l)  or  B(l)  if  a change  has  been  made. 

In  order  to  perform  arithmetic  with  the 
IU  and  IV  a L0GICAL*1  array  IW(- ) of  dimension  4 is  defined  and 
equivalenced  with  an  INTEGER*!  variable  IVAL.  IU  and  IV  are  thus 
placed  in,  or  extracted  from,  the  rightmost  byte  of  IVAL  by  passing 
them  through  IW (4 ) . 


Parameters  Passed 

MC12A  requires  as  parameters  the  arrays  A('),  IAv' ) and 
ISC(.,2)  which  are  in  blank  common  for  all  other  subroutines.  An 
array  IP(.,2^  with  first  dimension  NRMAX+1  is  also  passed.  This 
array  is  equivalenced  to  IE(').  IP(l,l)  contains  IA(NR0W+l)  for 
1=1,  NROW+1 , the  pointers  to  the  structural  columns  (.the  M matrix), 
while  IP(l,2)  contains  only  the  integers  1 through  NROW+1.  A 
REAL+4  array  W£PACE(’  ) of  dimension  4+NRMAX,  euuivalenoed  to  Y("  ■ 
is  also  passed.  IP(.,2)  and  WSPACF.(‘)  must  be  redimensioned  if 
NRMAX  is  changed.  , 


Method 


i 


, i 


l/  First  the  sum  of  squares  of  exponents  of  the  M matrix  is 
computed  and  printed. 

2/  The  array  IP(.,2)  is  set  up  and  MC12A  is  called. 

3/  If  MC12A  detects  an  empty  row  or  column  scaling  is  abandoned. 

b/  After  MC12A  has  returned  the  row  and  column  scales  for  M 
in  the  array  ISC (.,2)  the  structural  column  in  A(" ) and  the 
right-hand  side  3('  ) are  scaled  by  modifying  their  exponents. 

5/  The  new  sum  of  squares  of  the  exponents  of  the  one  matrix  is 
computed  and  printed. 

Messages 

MATRIX  SCALING  CALLED  FOR 

IFSCAL  was  set  equal  to  a nonzero  value  and  M will  be  scaled. 


INITIAL  SUM  0?  SQUARES  OF  EXPONENTS  = 'number' 

Initial  status  of  matrix  before  scaling  by  exponent. 

NEW  SUM  OF  SQUARES  OF  EXPONENTS  = ’number’. 

Status  of  matrix  after  scaling.  All  scaling  is  done  on  exponents 
only,  and  hence  scaling  factors  are  powers  of  lb.  If  the  matrix  was  well 
scaled  to  begin  with,  the  new  sum  of  squares  may  be  larger  than  the  old. 


i; 

fl 


•• 

i 


I 


I 


SINGULAR  M MATRIX-NO  SCALING  PERFORMED 


The  scaling  routine  detected  a zero  row  or  column.  The  run  will 
continue  but  normal  termination  is  unlikely. 

SUBROUTINE  CALLED 

MC12A 


SHFTE 

Shifts  the  etas  corresponding  to  the  upper  triangular  factor  (u) 
produced  by  INVERT  when  it  carries  out  an  LU  decomposition,  so  that 
they  are  contiguous  with  the  etas  for  the  L factor.  The  unified  eta 
file  is  then  ready  for  standard  product  form  application  and  update. 

Formal  Declaration 

SUBROUTINE  SHFTE 

Method 

The  U (backward'1  etas  were  built  up  in  reverse  order  from  the 
bottom  of  the  LE1.-  IE (*  ) and  E(-  regions  starting  at  NTMAX+1, 

NEMAX  and  NEMAX,  respectively.  The  number  of  U etas  (NUETA)  and 
eta  elements  (NUELEM)  is  known.  Given  NLETA  and  NLELEM  the  U etas 
are  shifted  to  follow  on  directly  from  the  L etas. 

SHIFTR 

Shifts  the  contents  of  one  region  to  another  (.completely  over- 
writing the  first  NROW  elements'. 


Formal  Declaration 


SUBROUTINE  SHJ.FTR  (I0LD,  INEW ) 


Parameters:  I$LD  (input):  Region  number  to  be  copied  (shifted) 

INEW  (input):  Region  number  into  which  region  I$LD  is 

to  be  shifted  (copied;. 


Notes 


(. ),  YTEMP(. ) are  numbered  1,  2 


and  k for  the  purposes  of  this  routine  and  the  parameters 


IOLD,  INEW  must  take  values  1,  2,  3>  or  4 


2 l The  array  BARRAY ( • ),  equivalenced  to  B(’ ) must  be  given  the 


dimension  4*NRMAX  if  NRMAX  is  charged 


UNPACK 


Unpacks  a packed  vector  in  the  coefficient  matrix  (a)  into 


the  work  region  Y (*  ) 


Formal  Declaration 


SUBROUTINE  UNPACK  (IV) 


Parameters:  IV (input) — the  secuence  number  (column  number)  of 


the  column  to  be  unpacked 


Method 


The  vector  is  unpacked  exactly  as  described  in  the 


Matrix  Storage 


section 


•JPBETA 

Updates  the  solution  voclc-  ” * ' when  ' v.'p'  is  pivoted 

on. 


iormal  Deelarat  ion 

SUBROUTINE  UPSET  A 


Method 

Standard  pivoting  procedure 
l/  Compute  new  DP  - X(.TR0WP)/y(lP.OWP) 

2/  Update  X(l)  - X(l)  - V,I  )*DP  for  I = i,moW 

Replace  X(IROWP)  = DP 

WRETA 

Write  a new  eta  to  the  end  of  the  eta  " le  using  the  vector 
in  Y(’)  and  IROWP. 

Formal  Dec It  rati on 

SUBROUTINE  WRETA 


Method 


1/  Increase  NEIEM  by  1 and  write 
TE(NELEM)  = IROWP 
F (NEIEM ) = Y( IROWP' 


(see  section 


rans  format  ion  ! roces  sin. 


2/  Zero  out  YlIROW?) 

3/  Pick  up  all  elements  in  '{[’)  with  absolute  value  greater  than 
ZTETA  and  pack  sequentially  in  IE(*)  and  E('  , increasing 
NELEM  at  each  step 

4/  Increase  NETA  by  1 and  define  EE  (NETA+l ) = NEEEM+1 

XCOIS 

Constructs  extra  columns  to  make  up  a square  matrix  when  NQUAD 
has  been  specified  and  only  NQUAD  (<  NRQW)  columns  were  supplied. 

The  matrix  M is  assumed  to  correspond  to  a quadratic  program  and 
be  of  the  form 

T 

D -A 
A 0 

XCOLS  picks  up  the  last  NROW-NQUAD  rows  of  the  matrix  supplied  and 

r -at  i 

constructs  I I , which  columns  are  added  on  to  the  assumed 

j to  give  the  full  M.  No  empty  rows  of  A are  allowed. 

Formal  Declaration 

SUBROUTINE  XCOI£ 


Method 

Search  rows  NQUAJ>1  through  NROW  in  sequence  for  nonzeros. 
Store  their  negatives  as  new  columns. 


5? 


Messages 


NUMBER  OF  ELEMENT  EXCEEDS  IAMAX-AB  • 


The  new  elements  er.  ee d storage  c-paj:'y. 


TRANSPOSE  COLUMNS  ADFED  ?(  PIATT  X 


Successful  completion 


' ' SOW  'row  no. ' ENCO  ED-ABANDONING  PROBLEM 

An  empty  row  (number  'row  no.')  has  been  encountered  in 
T 

trying  to  construct  -A  , Empty  columns  are  not  allowed.  Fatal  error. 
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